%% First load data arrays for analysis
clear;
fnames = dir("monod*.mat");
for k=1:length(fnames)
    load(fnames(k).name);
end

%% Define molar mass
% Define approximate molar mass of components
molmass_casamino = 135.5; % g/L

% Total carbon molar mass
molmass_NZCasePlus = 125/molmass_casamino;
clc;
% close all;

%% Define fitting function
monodfit = @(b,x) b(1).*x./(b(2) + x);
x0 = [1 0.01];
T = [25 30 37]';

Km_total = [];
Km_total_err = [];

figure;
%% Do 25C
tempmat = monod_mat_NZCasePlus_25C;

c = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
g = nanmean(tempmat(:,2:4),2);
ge = nanstd(tempmat(:,2:4),0,2);

f = fitnlm(c,g,monodfit,x0, 'Weights', 1./ge.^2);

Km = f.Coefficients.Estimate(2);
Km_err = f.Coefficients.SE(2);

cplot = linspace(0,max(c),1000);
g_estimate = monodfit(f.Coefficients.Estimate,cplot);

Km_total = [Km_total; Km];
Km_total_err = [Km_total_err; Km_err];

% Do fit on corrected concentration (in mM)
errorbar(c,g,ge, 'bo', 'linewidth', 1);
hold on;
plot(cplot,g_estimate,'b','linewidth',1);


%% Do 30C
tempmat = monod_mat_NZCasePlus_30C;

c = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
g = nanmean(tempmat(:,2:4),2);
ge = nanstd(tempmat(:,2:4),0,2);

f = fitnlm(c,g,monodfit,x0, 'Weights', 1./ge.^2);

Km = f.Coefficients.Estimate(2);
Km_err = f.Coefficients.SE(2);

cplot = linspace(0,max(c),1000);
g_estimate = monodfit(f.Coefficients.Estimate,cplot);

Km_total = [Km_total; Km];
Km_total_err = [Km_total_err; Km_err];

% Do fit on corrected concentration (in mM)
errorbar(c,g,ge, 'ko', 'linewidth', 1);
hold on;
plot(cplot,g_estimate,'k','linewidth',1);

%% Do 37C
tempmat = monod_mat_NZCasePlus_37C;

c = tempmat(:,1)/100*molmass_NZCasePlus*1000; % Convert to mM
g = nanmean(tempmat(:,2:4),2);
ge = nanstd(tempmat(:,2:4),0,2);

f = fitnlm(c,g,monodfit,x0, 'Weights', 1./ge.^2);

Km = f.Coefficients.Estimate(2);
Km_err = f.Coefficients.SE(2);

cplot = linspace(0,max(c),1000);
g_estimate = monodfit(f.Coefficients.Estimate,cplot);

Km_total = [Km_total; Km];
Km_total_err = [Km_total_err; Km_err];

% Do fit on corrected concentration (in mM)
errorbar(c,g,ge, 'ro', 'linewidth', 1);
hold on;
plot(cplot,g_estimate,'r','linewidth',1);

%% Set figure properties
set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')
box off;

%% Plot Km vs. temperature
figure;
hAx=axes;
errorbar(T,Km_total, Km_total_err, 'ko-', 'linewidth',1, 'markerfacecolor','g');
hAx.YScale='log';
xlim([23 39]);
ylim([.02 10]);
yticks([.02 .1 .5 1 3 10])
set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('T (°C)')
ylabel('Km')
box off;
legend 'Casamino acids';

otherplots = 0;
if otherplots
    figure;
    errorbar(T,Km_total, Km_total_err, 'ko-', 'linewidth',1);
    xlim([23 39]);
    set(gcf, 'position', [0 0 400 300])
    set(gca, 'fontsize', 20);
    xlabel('Temperature (°C)')
    ylabel('K_M (mM)')
    box off;
    
    figure;
    errorbar(1./(T+273.15), log(Km_total), Km_total_err./Km_total, 'ko-', 'linewidth',1);
    xlim(1./[273+39 273+23]); 
    ylim([-4 3])
    set(gcf, 'position', [0 0 400 300])
    set(gca, 'fontsize', 20);
    xlabel('1000/T (K-1)')
    ylabel('Log(Km)')
    box off;
end